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We study the competition of magneto-dipole, anisotropy and exchange interactions in compos¬ 
ite three dimensional multiferroics. Using Monte Carlo simulations we show that magneto-dipole 
interaction does not suppress the ferromagnetic state caused by the interaction of the ferroelectric 
matrix and magnetic subsystem. However, the presence of magneto-dipole interaction influences the 
order-disorder transition: depending on the strength of magneto-dipole interaction the transition 
from the ferromagnetic to the superparamagnetic state is accompanied either by creation of vortices 
or domains of opposite magnetization. We show that the temperature hysteresis loop occurs due to 
non-monotonic behavior of exchange interaction versus temperature. The origin of this hysteresis is 
related to the presence of stable magnetic domains which are robust against thermal fluctuations. 
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I. INTRODUCTION 

Multiferroic materials are materials with coupled mag¬ 
netic and electric degrees of freedom. 1 5 One example of 
this coupling is due to spin-orbit interaction in certain 
crystals. However, this coupling is relatively weak. Cur¬ 
rently there is an active search for multiferroic materials 
with strong couplingP^ 

One way to strongly coupled magnetic and electric 
degrees of freedom is to develop hybrid ferroelectric- 
ferromagnetic layered materials where mechanical stress 
produces strong correlations between the layersPQ 1 Re¬ 
cently another promising possibility has been suggested 
based on granular materials where small metallic ferro¬ 
magnetic (FM) grains were embedded into ferroelectric 
(FE) matrix or these grains were located in close proxim¬ 
ity to the FE substrate.^ The presence of small metallic 
grains increases the strength of Coulomb interaction pro¬ 
viding the necessary coupling between the FE and FM 
degrees of freedom. 

An important question is to understand the nature 
of multiferroic state in granular multiferroic materials. 
On the mean-field level the properties of these materials 
have been understood.^ It was shown that the exchange 
coupling J depends on the properties of the FE ma¬ 
trix/substrate, in particular on the dielectric permittiv¬ 
ity e of the surrounding medium. 12 ^ Due to temperature 
dependence of dielectric permittivity, e(T), the exchange 
interaction depends non-monotonically on temperature 
leading to the inverse phase transition with paramag¬ 
netic phase appearing at lower temperatures compared 
to the ferromagnetic phase. The qualitative behavior of 
J is different for small and large (compared to the grain 
size 7g r ) inter-grain distances a: for large inter-grain dis¬ 
tances, a > r gr , the J-value is increased in the vicinity of 
the FE Curie point due to suppression of the Coulomb 


blockade effects leading to a different magnetic state at 
these temperatures. 

However, it is still an open question whether the 
magneto-electric coupling obtained in the mean-field the¬ 
ory is robust against the magneto-dip ole and anisotropy 
interactions neglected in the mean-field approach. We 
use the numerical modelling to address this question. 

We investigate the non-equilibrium meta-stable states 
in granular system and study the nature of new rneta- 
stable phases appearing in the system with temperature 
dependent exchange interaction. In addition, we answer 
the question if temperature dependent exchange interac¬ 
tion can lead to unusual blocking effects. 

To be more specific, we study the magnetic behavior 
of composite multiferroics with T™ < T™, where T™ 
is the FE Curie temperature of paraelectric-ferroelectric 
transition of the FE matrix and T™ is the Curie tem¬ 
perature of FM grains. We focus on the temperature 
range T <C T™ where all grains are in the FM state 
and study the phase diagram of granular multiferroics 
beyond the mean-field approximation using Monte-Carlo 
simulations. In particular, we study the combined effect 
of magnetic anisotropy, the long-range magneto-dipole 
(MD) interaction and the exchange coupling. We show 
that there is an inverse magnetic transition in the sys¬ 
tem which is robust with respect to the MD interaction. 
This inverse transition disappears only for strong MD in¬ 
teraction, stronger than the exchange coupling. In addi¬ 
tion, we show that non-monotonic temperature behavior 
of intergrain exchange interaction leads to a new type of 
hysteresis in composite multiferroics. 

Three-dimensional (3D) nanostructures composed of 
single-domain ferromagnetic particles has been in¬ 
tensively studied both experimentally and theoreti¬ 
cally.' 1422 The interplay of magnetic anisotropy, long- 
range magneto-dipole interaction and short range ex- 
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change interaction defines the magnetic state of the sys¬ 
tem. Depending on the ratio of these interactions differ¬ 
ent magn etic states are possible in granular ferromag- 
netsPSH) Among them are superparamagnetic (SPM), 
super spin-glass (SSG) and superferromagnetic (SFM) 
states. 

The most studied situation is related to the case of 
large intergrain distances (> 2nm) and small exchange 
interaction where magnetic state is defi ned by t he compe¬ 
tition of MD interaction and anisotropy! 21 * 25 * 2 ^ * The mag¬ 
netic anisotropy is responsible for “blocki ng” p henomena 
and defines the blocking temperature Tb.^ 8 ^* The weak 
MD interaction modifies the blocking temperature, while 
the strong MD interaction leads to the SSG state“H2H 

For small distances between grains (~ 1 nm) the ex¬ 
change interaction is c rucial. It leads to the formation 
of the SFM statePlMM! In such system s the SPM-SFM 
transition occurs. Even a weak exchange interaction can 
influence the magnetic state of the system.® In particu¬ 
lar, the FM ordering with long range ferromagnetic cor¬ 
relations appears. 

Different theoretical methods have been used to study 
granular magnets. The mean-field approach allows to 
study granular systems with finite short range exchange 
interaction and zero (or weak) MD interaction neglect¬ 
ing fluctuation effectsW® Modelling based on Landau- 
Lifshitz equation allows to consider the MD, magnetic 
anisotropy and exchange interactions.® However, this 
approach needs to be generalized in the presence of ther¬ 
mal fluctuations by introducing the Langevin forces.® 
These fluctuations are important for granular magnetic 
systems since the granular magnetic moment is rela¬ 
tively small and fluctuation effects are pronounced es¬ 
pecially near the phase transitions. However, the inclu¬ 
sion of Langevin forces in the nonlinear spin dynamics 
is not numerically efficient. Therefore we use Monte- 
Carlo (MC) simulations which allow to study phase 
transitions in composite multiferroics with strong long- 
range MD interaction and arbitrary thermal fluctua¬ 
tions ! 141151211241261411 ® 


MC modelling strongly depends on the degree of 
anisotropy. At strong anisotropy the problem reduces 
to the Ising model with magnetic moment of each parti¬ 
cle having only two directions defined by the anisotropy 
axis. In this case the MC modelling is very efficient and is 
based on the trial spin flips. Another type of MC mod¬ 
elling is based on the Heisenberg model with arbitrary 
magnetization direction. This model is more general but 
it requires more computational time due to sp in rota¬ 
tions over the sphere rather than spin flips IA2IL 4 I 46 I 47 I We 
use our own MC code with random spin-flips and random 
spin-rotations that is valid for any anisotropy. 


The paper is organized as follows: In Sec. [IT] we formu¬ 
late our main results. In Sec. m we discuss the model of 
composite multiferroics. In Sec. IV we introduce impor¬ 
tant physical quantities which we calculate. We discuss 
our results in Sec. El The details of our numerical calcu¬ 
lations are presented in Appendix [A] 


II. MAIN RESULTS 


Here we summarize our main results. The non¬ 
monotonic temperature dependence of exchange inter¬ 
action in composite multiferroics leads to the unusual 
evolution of the magnetic state with temperature. The 
intergrain exchange interaction has either peak or deep 
in the vicinity of the FE phase transition due to cou¬ 
pling of electric and magnetic degrees of freedom. In the 
mean field approximation the peak in the exchange in¬ 
teraction leads to the onset of FM state in the vicinity 
of FE phase transition. The deep in the exchange in¬ 
teraction suppresses the FM state in the vicinity of the 
FE Curie point. We use Monte-Carlo simulations to in¬ 
vestigate the influence of long-range MD interaction and 
magnetic anisotropy on the magnetic phase diagram of 
composite multiferroics. We show that MD interaction 
and anisotropy do not suppress the magneto-electric cou¬ 
pling in these materials, however their interplay produces 
a new type of hysteresis. Our results are the following: 

1) The Monte-Carlo simulations reproduce the mean 
field results in the absence of MD interaction and mag¬ 
netic anisotropy. Similar to the mean field approach, the 
FM state exists in the vicinity of the FE Curie point and 
the disordered state appears away from this region. 

2) The finite MD interaction does not suppress the FM 
ordering in the vicinity of the FE phase transition even if 
the MD interaction is twice stronger than the exchange 
interaction. The presence of MD interaction leads to the 
appearance of the domain structure and to splitting of 
uniform FM state. This result is similar to Ref. [34, where 
a weak FM interaction leads to the formation of FM do¬ 
mains, while strong MD interaction produces a vortex 
structure. 

3) The magnetic state depends on the strength of MD 
interaction outside the FM region: the system is in the 
SPM state for weak MD interaction and in the antiferro¬ 
magnetic stripe phase for strong MD interaction. 

4) The magnetic anisotropy does not influence the FM 
state. However, it prevents the formation of vortices in 
the transition region and leads to a widening of the FM 
domain. 

5) The “blocking phenomenon” does not appear at 
finite magnetic anisotropy and zero MD interaction at 
considered temperatures meaning that the system has 
enough time to reach the ground state such that the non¬ 
equilibrium state does not appear. 

6) The “blocking phenomenon” appears at finite mag¬ 
netic anisotropy and finite MD interaction. The temper¬ 
ature hysteresis loop occurs due to non-monotonic be¬ 
havior of exchange interaction vs. temperature. The 
origin of this hysteresis is related to the presence of sta¬ 
ble magnetic domains which are robust against thermal 
fluctuations, MD, exchange, and anisotropy interactions. 

7) The AFM stripes appear in the case of deep in the 
exchange interaction. 

Below we discuss these results in details. 
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III. THE MODEL 
A. Magnetic subsystem 

In this subsection we consider the magnetic subsys¬ 
tem. We model a composite multiferroic as an ensemble 
of FM grains embedded into FE matrix. All grains are 
homogeneously magnetized single domain FM particles 
of the same volume V and saturation magnetization M s . 
For temperatures T -C T™ , the saturation magnetiza¬ 
tion M s is constant. Each grain with volume V has a 
magnetic moment g = M S V and is treated as a point 
dipole located at the centre of the grain. The grains are 
pinned to the sites of the regular cubic lattice with lattice 
spacing a and can freely rotate their adjusting magnetic 
moments. 

The whole system is modelled as a 3 D lattice of clas¬ 
sical spins, with magnetic moment of the ith grain being 
Hi = /iSwhere the unit vector S* = (Sf,Sf,Sf) is 
the spin of ?’th particle representing the direction of the 
magnetic moment. 

We assume that each grain has a uniaxial anisotropy. 
Spatial distributions of anisotropy axes varies in different 
experiments and depends on the preparation condition. 
The anisotropy axes can be homogeneously distributed 
over the solid angle, or uniformly distributed in a certain 
plane. We assume that the easy axes of all grains are 
oriented in the z-direction. This situation is realized in 
experiment with magnetic field applied during the sample 
preparation.^ 

The Hamiltonian of the system has the form 

T~i = Aexc + Adip + 'Ha.n- (1) 

The first term, "H exc , describes the exchange coupling 
between grains i and j 

U exc = ~jJ2^- S v ( 2 ) 

(i,j) 

where the sum is over the nearest neighbour pairs of 
grains. 

The second term, in Eq. 0 describes the long- 

range magneto-dipole (MD) interaction between mag¬ 
netic moments Hi and l^j °f individual grains 

t, • Sj rf ■ — 3(Si • Tij)(Sj ■ v.ij) 

^dip =9 2 ^ -Is-> (3) 

i<j b 

where is the distance between magnetic moments at 

sites i and j measured in units of lattice spacing a and g 
is the MD interaction constant. 

The third term, "H an , in Eq. 0 describes uniaxial 
anisotropy energy 

H an = -A'^(e„,-S l ) 2 , (4) 


where K is the temperature independent magnetic 
anisotropy energy of a single grain. The unit vector e~ 
defines the direction of the anisotropy easy axis. 

We consider the energy parameters (J, g , K, T) in arbi¬ 
trary units. Parameters, g and K depend on a grain vol¬ 
ume and can be controlled by varying the grain size. The 
dipole coupling g is additionally depends on the lattice 
spacing a. This allows one to vary parameters g and K in 
a wide range. The exchange interaction is proportional 
to the grain surface and scales with the volume as R 2 / 3 . 
Moreover, the ratio of MD and exchange interactions can 
be controlled by varying the interparticle distance a. The 
MD interaction decays as 1/a 3 with distance, while ex¬ 
change interaction decays exponentially e ~ Ka , where n is 
the inverse length which depend on the band structure 
of the surrounding FE matrix and the Fermi energy of 
electrons inside grains. 

B. Ferroelectric subsystem 

The FE matrix is characterized by the Curie tempera¬ 
ture T EE . We study the temperature region in the vicin¬ 
ity of T EE . The most important characteristic of the 
FE matrix in our consideration is the FE dielectric per¬ 
mittivity e which has a peculiarity in the vicinity of the 
phase transition point, T EE PSlThis peculiar behavior of e 
provides a strong co uplin g between electric and magnetic 
degrees of freedom . 12 13 

C. Interaction of magnetic and ferroelectric 
subsystems 

Recently the intergrain exchange interaction, J , in 
composite multiferroics was s tudied and its temperature 
behavior was predicted . 12 ^ In the vicinity of the FE 
phase transition the exchange interaction has a peculiar¬ 
ity. Depending on the system parameters it has either 
peak or deep. Such peculiarity appears due to combine 
influence of Coulomb blockade effects and the tempera¬ 
ture dependence of dielectric permittivity of the FE ma¬ 
trix on the exchange interaction. The peculiarity of the 
exchange interaction in the vicinity of the FE Curie point 
T = T EE is related to the peculiarity in the dielectric per¬ 
mittivity e of the FE matrix. The exchange interaction 
as a function of e has the form 12 

J(T) = J 0 e^ a/r ^-\ (5) 

where parameter Jo > 0 is the e-independent part of 
exchange interaction, Jo decays exponentially with inter- 
grain distance a; 7 is the numerical coefficient of order 
one. The dielectric permittivity e has a peak at tempera¬ 
ture T = T EE . For 7 < r gr /a (small intergrain distances) 
the exchange interaction has a deep, while for 7 > r gr /a 
(large intergrain distances) it has a peak. 

The dependence of the intergrain exchange interaction 
J on the FE permittivity e is the signature of magneto- 
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electric coupling emerging in composite multiferroics. 
The peak of exchange interaction leads to the unusual 
magnetic phase diagram: the FM state appears in the 
vicinity of the FE Curie temperature T EE , while away 
from T EE the system is in the SPM state. Thus, the FM 
state in the system exists in a finite region around FE 
Curie point, T EE . In case of deep in the temperature de¬ 
pendence of exchange constant J the opposite situation 
occurs: the FM state is suppressed in the vicinity of T EE 
due to interaction of magnetic and FE subsystems. 

Above effects have been studied in Ref. □J using the 
mean field approximation without taking into account 
the MD interaction and anisotropy. Here we take into 
account both interactions and study the influence of MD 
interaction and magnetic anisotropy on the magneto¬ 
electric coupling in composite multiferroics. 

For large inter-grain distances, 7 > r gr /a, the exchange 
interaction J(T ) has the peak. In this work we model this 
peak as follows 

J(T) = J 0 e-( T - T c E ) 2 /- 2 , ( 6 ) 

where w is the width of the exchange peak and Jo is 
the amplitude of intergrain exchange interaction. The 
peak in J(T) occurs at temperatures T = T EE , when the 
Coulomb blockade is suppressed and the electron wave 
functions become weakly localized leading to a strong 
overlap and strong exchange interaction.^ The Coulomb 
blockade is suppressed when permittivity e(T) reaches its 
maximum value at T = T EE . 

For small inter-grain distances, 7 < r gr /a, the ex¬ 
change interaction J(T) has a deep which we model as 
follows 

J(T) = J 0 (l-e- (T - T £ E)2/u;2 ) . (7) 


where summation is over the all grains in a given cell 
N c = L j? and the averaging is defined as summation over 
the all possible positions of the cell centre, () = N~ x JA. 

Next, we introduce the spin-spin correlation function G 
as an averaged correlation function for nearest-neighbour 
pairs 

G = W = sE (s -' s j!. (io) 

R g (i,j) 

where g are the six vectors of nearest neighbours and the 
last summation (i, j) is over the all nearest-neighbour 
pairs in the lattice. The correlation function G is im¬ 
portant for understanding magneto-transport in granu¬ 
lar magnets. The magneto-resistance (MR) of granular 
magnetic film is proportional to this correlation function, 
MR(T)~ G. The MR measurements can be considered 
as the probes of the magnetic state of the system. 

V. DISCUSSION OF RESULTS 

A. Influence of magneto-dipole (MD) interaction 
on the magnetic phase diagram of composite 
multiferroics 

In this subsection we discuss the influence of MD in¬ 
teraction on the magnetic phase diagram of composite 
multiferroics for the case of large inter-grain distances, 
where exchange interaction J has a peak around T EE , 
see Eq. (| 6 |. 

1. Zero magneto-dipole (MD) interaction 


IV. CALCULATED QUANTITIES 


In this section we introduce physical quantities which 
we calculate. The calculation procedure is described in 
Appendix [A] The first quantity we calculate is the aver¬ 
age magnetizatiorP^^ 


M(T) = 


1 

N 



( 8 ) 


where N is the number of lattice sites. 

In the presence of MD interaction and zero external 
magnetic field, the average magnetization M(T) is not 
an efficient quantity: the lattice spins form complex mag¬ 
netic patterns, either domains or vortices, and the mean 
magnetization vanishes, 1 /IVjASi = 0 , even if locally 
the system is in the FM state. To account for local FM 
correlations we introduce the cell averaged magnetiza¬ 
tion, m(T), over a cell with linear size L c . We find that 
the optimal size for such cell average is L c = 5 


m(T) = 



(9) 


In the absence of long-range MD interaction and 
anisotropy the magnetic phase diagram obtained using 
the Monte-Carlo simulations coincides with the mean- 
field phase diagram (see Fig. [I]). We used the following 
parameters: J 0 = 2.5, AJ FE = 0.1, T EE = 0.5. These 
parameters correspond to Fe grains embedded into or¬ 
ganic TTF-CA FE matrix. The grain size is 5—10 nm and 
the intergrain distance > 1 nm. For these parameters the 
intergrain exchange interaction is about J ~ 300K. The 
Curie temperature of bulk ferroelectric TTF-CA is 80 K. 
However, for granular array it can be smaller; the Curie 
temperature of TTF-CA in composite granular metal/FE 
system is 50 K.^ 

The FM state appears around T EE , where J(T) has a 


3 ) and the SPM state 
^,) . The finite FM re- 
lase transitions due to 


maximum (see snapshot in Fig. D 
exists outside this region (see Fig. 
gion appears with two magnetic p 
peak in the exchange interaction. In the mean-field ap¬ 
proximation both transition temperatures T ,(f are defined 
as T = 2 J(T). Due to zero MD interaction domains are 
not formed, the ground state of the system is the homoge¬ 
neous FM state and the magnetization M ( T ) is coincide 
with the cell averaged magnetization m(T). In the SPM 
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FIG. 1. (Color online) Magnetic phase diagram of compos¬ 
ite multiferroic vs. temperature for zero magneto-dipole in¬ 
teraction (g = 0) and magnetic anisotropy (K = 0). Solid 
(red) line shows the temperature dependence of intergrain ex¬ 
change interaction, J{T). Straight dash dotted (orange) line 
stands for temperature T. Gray dash dot-dotted and green 
dotted lines show the average magnetization M(T ) and cell 
averaged magnetization m(T), respectively. Blue dashed line 
shows the nearest neighbour correlation function G. Transi¬ 
tion temperatures are defined using the mean-field equa¬ 
tion, T = 2 J(T). 



FIG. 2. (Color online) Snapshots of single magnetic layer 
of composite multiferroic. Panels (a) and (b) show the dis¬ 
ordered (SPM) magnetic state and the ordered FM state, re¬ 
spectively. Position of these panels is shown in Fig. 1 by black 
arrows. 

state the saturation magnetization and the correlator G 
tends to zero since the system is in the disordered state 
due to thermal fluctuations. 


2. Weak and moderate magneto-dipole (MD) interaction 

The long-range MD interaction competes with inter¬ 
grain exchange interaction suppressing the FM state in 
the system. Figure [3] shows the case of weak MD in¬ 
teraction with the dipole constant, g = 0.5. This value 
of g is typical for Fe grains with size a = 4 nm, where 



FIG. 3. (Color online) Magnetic phase diagram of compos¬ 
ite multiferroic vs. temperature T for weak magneto-dipole 
interaction (g = 0.5) and zero magnetic anisotropy (K = 0). 
All notations are defined in Fig. 1. 



FIG. 4. (Color online) Snapshots of single magnetic layer of 
composite multiferroic. Panel (a) shows the FM state divided 
into stripe domains appearing due to the interplay of exchange 
and weak MD interactions. Panel (b) shows the vortex state 
with magnetic vortices appearing due to strong MD interac¬ 
tion. Panel (c) shows FM chains forming antiferromagnetic 
pattern appearing at low temperatures and MD interaction 
being stronger than exchange interaction. Position of panels 
(a) and (c) is shown in Fig. 3 by black arrows. 


g = (2.2 /i b V/X F e )V a3 ~ 50 K, Ap e = 0.28 nm is the Fe 
lattice parameter.^ For these parameters the FM region 
exists. However, the MD interaction reduces the size of 
FM state and leads to the formation of domains in the 
system (see the left panel in Fig. [4]). Above the FM region 
the SPM state appears, similar to the case of zero MD in¬ 
teraction, meaning that the thermal fluctuations exceed 
the exchange and MD interactions. Below the FM region 
the thermal fluctuations exceed the exchange interaction 
but not the MD interaction. As a result the antiferromag¬ 
netic stripes appear at low temperatures, T < 0.35 (see 
panel (c) in Fig. [4|. The transition to stripe structure 
occurs via formation of large antiferromagnetic domains 
with temperature dependent sizes. 
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FIG. 5. (Color online) Magnetic phase diagram of compos¬ 
ite multiferroic vs. temperature T for intermediate magneto¬ 
dipole interaction g = 2.0 and zero magnetic anisotropy 
(K = 0). All notations are defined in Fig. 1. 


3. Strong magneto-dipole (MD) interaction 


For strong MD interaction, g > 2.0, the uniform FM 
state does not appear in the system (see Fig. [5]). How¬ 
ever, the long-range magnetic correlations are still exist 
in the system close to T EE . For such a strong MD interac¬ 
tion the vortex-like structure appears in the system (see 
the central snapshots in Fig.[4|b)). The vortex structure 
transforms into the stripe structure outside the T EE re¬ 
gion. For strong MD interaction the SPM state does not 
appear for temperatures above and below T EE , instead 
the stripe structure appears in these regions. At higher 
temperatures the stripe structure transforms into SPM 
state due to thermal fluctuations. 

The MD interaction grows with particle volume as V 2 . 
For Fe grains of size 6nm and interparticle distance lnm 
the dipole constant is g ss 230 K. This value of MD in¬ 
teraction equals to the peak value of exchange interac¬ 
tion. The exchange interaction grows with grain surface 
as F 2 / 3 and it is intergrain distance dependent. 

Figure [6] shows the case of strong MD interaction with 
dipole constant g = 5 being twice larger than the peak 
value of the exchange interaction. This case is typical for 
lOnm Fe grains with MD constant g exceeding the room 
temperature. Even in this case the FM domains exist in 
the vicinity of the FE Curie point. The FM state in the 
vicinity of the FE phase transition is robust against the 
MD interaction leading to the unusual magnetic phase 
diagram of composite multiferroics due to ME coupling 
of ferroelectric and magnetic subsystems. 



FIG. 6. (Color online) Magnetic phase diagram of composite 
multiferroic vs. temperature T for strong magneto-dipole in¬ 
teraction, g = 5.0, and zero magnetic anisotropy, K = 0. All 
notations are defined in Fig. 1. 


B. Influence of magnetic anisotropy on the 
magnetic phase diagram of composite multiferroics 

In this subsection we discuss the influence of magnetic 
anisotropy on the magnetic phase diagram of composite 
multiferroics. The magnetic anisotropy in granular ma¬ 
terials is much stronger than in bulk magnets due to sur¬ 
face and shape anisotropy! 16 * 18 ! It plays an important role 
in formation of magnetic state of granular magnets. The 
magnetic relaxation time in the system of non-interacting 
particles exponentially depends on the ratio of anisotropy 
energy and temperature, r r ~ exp(A'/(fcBT)). At low 
temperatures the relaxation time becomes larger than 
the characteristic measurement time. At these tempera¬ 
tures the measured magnetic properties are the proper¬ 
ties of non-equilibrium or “blocked” state. The tempera¬ 
ture hysteresis of magnetic properties is the signature of 
“blocking” phenomenon. The Monte-Carlo (MC) calcu¬ 
lations in some way are similar to real experiment: sim¬ 
ulations start with a certain non-equilibrium state and 
the system “relaxes” to the equilibrium state via discrete 
steps during the simulations. If number of MC steps Nmc 
(which can be associated with measurement time if the 
attempt frequency of the system is known) exceeds a cer¬ 
tain value N r (which can be associated with relaxation 
time r r ) then the system relaxes to the equilibrium state 
during simulations. In the opposite case, the system is 
locked into some non-equilibrium magnetic state. 

Figure [7] shows the magnetic behavior of granular 
multiferroics with strong anisotropy, I\ = 6.0. This 
anisotropy is twice larger than the peak value of the ex¬ 
change interaction. This situation is realized for 6 nm Fe 
grains with anisotropy constant I\ = 0.8T0” 6 erg/cm 3 £20 
Figure [7] coincides with mean-field theory, where FM 
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FIG. 7. (Color online) Magnetic phase diagram of composite 
multiferroic vs. temperature T for zero magneto-dipole in¬ 
teraction and strong magnetic anisotropy ( I\ = 6.0). Insert: 
Snapshot of a single magnetic layer of composite multifer¬ 
roic with disordered magnetic state due to strong anisotropy. 
Black arrow indicates the position of this snapshot on phase 
diagram. All notations are defined in Fig. 1. 


state exists in the vicinity of the FE phase transition 
and the disordered state appears away from T EE . The 
only difference with mean field theory is related to the 
fact that magnetic moments in the disordered state have 
only two possible directions along the anisotropy axis (see 
inset in Fig. [7|). 

For zero MD interaction the MC results do not de¬ 
pend on the initial state of the system. We use two dif¬ 
ferent approaches for MC calculations, see the details in 
the Appendix: 1) The starting configuration is the FM 
alignment along the z-direction at the initial temperature 
point. After a certain number of MC steps the resulting 
spin configuration is used for the next temperature point. 
2) The starting configuration is the disordered state at 
each temperature point. In addition, we change the di¬ 
rection of temperature evolution and the number of MC 
steps. Both approaches lead to the same results without 
hysteresis behavior. Thus, we conclude that in our mod¬ 
elling Amc N t . And the magnetic anisotropy alone 
does not qualitatively change the magnetic phase dia¬ 
gram of composite multiferroics and it does not lead to 
the suppression of ME effects in the system. 


C. Hysteresis behavior of composite multiferroics 
with strong magnetic anisotropy and magneto-dipole 
interaction 

In this subsection we discuss the influence of strong 
magnetic anisotropy and MD interaction on the phase 
diagram of composite multiferroics and show that these 




FIG. 8. (Color online) Temperature hysteresis of composite 
multiferroic. Left and right panels correspond to the case of 
increasing and decreasing temperature, respectively. 


interactions lead to new features. We use the following 
approach: 1) We move from high to low temperatures 
starting from the FM state at initial temperature T = 
0.72. At each next temperature point we begin with the 
final state of the previous temperature point. 2) We move 
from low T = 0.28 to high T = 0.72 temperatures with 
the initial state corresponding to the FM state. Results 
are shown in Fig. [8j 

The presence of MD interaction increases the blocking 
temperature leading to larger magnetic relaxation time, 
7 ~ r ! 5 12 o | However, this is not the case for our system. To 
understand the nature of hysteresis in composite mul¬ 
tiferroics consider the starting point T = 0.28 in the 
“warming” case, left panel, with the uniform FM state 
as the initial state. The final state at this temperature 
obtained after MC simulations is the disordered (SPM) 
state. The system relaxes from FM to SPM state in 
the process of calculations. Therefore neither magnetic 
anisotropy nor the MD interaction produce “blocking” in 
the system at these temperatures. 

However, if we move from high to low temperatures 
the relaxation to the SPM state does not occur at tem¬ 
perature T = 0.28. The reason for the cooling process is 
related to the fact that one needs to pass the peak in the 
exchange interaction at temperature T = 0.5 where mul¬ 
tiple inhomogeneous FM state with domains of different 
orientations is formed. Around T = 0.5 the average mag¬ 
netization (the red curve) is zero. But the cell averaged 
magnetization is finite meaning that the system is divided 
into domains of opposite magnetization. These domains 
occur due to the interplay of exchange and MD interac¬ 
tions. This multidomain state is robust against thermal 
fluctuations for T < 0.5 even for zero exchange interac¬ 
tion. This is in contrast to the uniform FM state, which 
can be destroyed in the absence of exchange interaction 
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FIG. 9. (Color online) Evolution of magnetic phase diagram 
of composite multiferroic as a function of exchange interaction 
at fixed temperature. 


even at low temperatures. As a result the magnetization 
vs. temperature has a hysteresis behavior. 

At higher temperature, T ss 0.7, the multidomain state 
is unstable against thermal fluctuations at these temper¬ 
atures and the hysteresis behavior is absent. 

Figure [9] shows the magnetic phase diagram as a func¬ 
tion of exchange constant J at fixed temperature T = 0.3. 
On the left panel we increase the exchange constant J 
starting from zero to a certain high value and decrease it 
back to zero. The initial spin configuration at J = 0 is the 
uniform FM state. On the right panel in Fig. [9] we start 
with finite value of exchange interaction J, decrease it to 
zero and return back to the same value. Obviously, such 
a nonmonotonic behavior of exchange interaction occurs 
with changing the temperature in composite multiferroic. 
Figure [9] shows the hysteresis behavior caused by the non 
monotonic change of the intergrain exchange interaction. 
Similar hysteresis occurs as a function of temperature. 

To summarize, the above hysteresis is specific to com¬ 
posite multiferroics - materials with non-monotonic be¬ 
havior of exchange interaction. The hysteresis is absent 
for systems with temperature independent exchange in¬ 
teraction. The peculiar feature of this hysteresis is re¬ 
lated to the fact that it appears in the vicinity of the FE 
Curie point. 


D. Deep in the exchange interaction in the vicinity 
of FE Curie point 

The intergrain exchange interaction has either peak or 
deep in the vicinity of FE Curie point depending on the 
system parameters? 12 Here we study the situation with 
deep in the exchange interaction in the vicinity of the FE 


Curie point. Figure [TU^, shows the case of zero MD inter¬ 
action and zero anisotropy with the following parameters: 
the deep is J 0 = 0.75, A T FB = 0.07, and T™ = 0.5. In 
this case the Monte-Carlo simulations and the mean field 
approximation coincide. For chosen parameters the ex¬ 
change interaction exceeds the temperature in the whole 
range except the close vicinity of FE transition, where 
exchange interaction is small and the system is in the 
SPM state. Outside this region the system is in the FM 
state. A single domain state with two magnetic phase 
transitions is realized for zero MD interaction. 

A moderate MD interaction leads to the domains for¬ 
mation in the FM phase and to the formation of magnetic 
vortices in the transition regions. Figure |T7T} > shows this 
behavior for dipole constant g = 0.5. 

Strong MD interaction (g = 5.0) leads to the suppres¬ 
sion of the FM state at high temperatures, see Fig. pTi] '. 
Here the AFM stripe structure appears instead of FM or¬ 
dering. At low temperatures the FM state exists. There¬ 
fore in the case of deep in the exchange interaction the 
MD interaction leads to the suppression of high temper¬ 
ature magnetic phase transition in contrast to the case 
of peak in the exchange interaction. At higher tempera¬ 
tures, T > 1, the AFM state is suppressed due to thermal 
fluctuations. 


E. Applicability of results 

Here we discuss the applicability of our approach. 
First, we study the case of regular magnetic array with 
fixed intergrain distances. In real materials this distance 
fluctuates leading to the dispersion of MD and exchange 
interactions. 

Second, we consider 3D multiferroic materials which 
produced via bottom-up method. A different top-down 
fabrication, based on layer by layer growth, is used to 
produce a single layer of magnetic grains. In 2D systems 
the influence of MD interaction on the magnetic phase 
diagram is different from 3D case. This situation requires 
further investigation. 


VI. CONCLUSION 

We studied the competition of magneto-dipole, 
anisotropy and exchange interactions in composite three 
dimensional multiferroics - materials with magnetic 
grains embedded into FE matrix. The peculiarity of com¬ 
posite (or granular) multiferroics is related to the fact 
that interparticle interaction is affected by the FE ma¬ 
trix. Granular multiferroics show the magneto-electric 
coupling effect. Using Monte Carlo simulations we 
showed that magneto-dipole interaction does not sup¬ 
press the ferromagnetic state caused by the interaction of 
the ferroelectric matrix and magnetic subsystem. Thus, 
MD interaction does not suppress the ME effect in gran¬ 
ular multiferroics. However, the presence of magneto- 
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FIG. 10. (Color online) Magnetic phase diagram of composite multiferroic vs. temperature T for the case of deep in the 
exchange interaction in the vicinity of the FE Curie point. Panel (a): zero magneto-dipole (MD) interaction (g = 0.0) and 
zero magnetic anisotropy (K = 0.0). Panel (b): moderate MD interaction (g = 1.0) and zero magnetic anisotropy ( I\ = 0.0). 
Panel (c): strong MD interaction (g = 1.0) and zero magnetic anisotropy (K = 0.0). Solid (red) line shows the temperature 
dependence of intergrain exchange interaction, J(T). Straight dash dotted (orange) line stands for temperature T. Gray dash 
dotted and green dotted lines show the average magnetization M(T) and cell averaged magnetization m(T), respectively. Blue 
dashed line shows the nearest neighbour correlation function G. 


dipole interaction influences the order-disorder transi¬ 
tion: depending on the strength of magneto-dipole in¬ 
teraction the transition from the FM to the SPM state 
is accompanied either by creation of vortices or domains 
of opposite magnetization. 

We showed that “blocking phenomenon” appears at fi¬ 
nite magnetic anisotropy and finite MD interaction. The 
temperature hysteresis loop occurs due to non-monotonic 
behavior of exchange interaction vs. temperature. The 
origin of this hysteresis is related to the presence of sta¬ 
ble magnetic domains which are robust against thermal 
fluctuations. 
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Appendix A: Calculation procedure 

We use classical Monte Carlo (MC) simulations and 
the standard Metropolis algo rithm to model magnetic 
properties of the system E3E2H52] We consider L x L x L 
(L = 20) cubic lattice with periodic boundary conditions. 
To efficiently evaluate the long-range MD interaction in 
systems with relatively small number o f part icles (as, for 
example, L = 5, 6, 7 considered in Ref. 56 60 ) one has to 
implement Ewald summation technique! 61 * 62 * We account 
the MD interaction by direct summation in the real space 
applying the minimum image convention. 6 ^ In terms of 
the range of the interaction that have been taken into ac¬ 
count, this scheme is equivalent to the fast Fourier trans¬ 
form method used for micromagnetic simulations ! 63 * 64 ^ 

We use the FM state ordered along the ^-direction, 
Si = b z , as an initial spin configuration for simulating 
at the first temperature point. The resulting spin state 
is used as an initial state for the next temperature point 
and so on. To study hysteresis effects we make two pas¬ 
sages: first, we start with low temperature and increase 
the temperature during our calculations; second, we do 
the opposite. 

One MC step consists of Lx Lx L consecutive changes 
in the lattice spin orientations. We calculated the change 
in the energy of the system AH'. if change is negative, 
A H < 0, a new state is accepted; if change is posi¬ 
tive, AH > 0, the new state is accepted with probability 
e-A -H'T. In our simulations we use Nmc = 12000 MC 
steps per spin and study 60 samples in every 200 MC 
steps to calculate thermal properties. To check the sta¬ 
bility of final configuration on the number of MC steps 
we increased the number of MC steps five times (up to 
Nmc = 6 ■ 10 4 ) and find no difference in the resulting 
state. 
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We generate the update in spin directions using two 
ways. First, the spin orientations were distributed uni¬ 
formly over the unit sphere’s surface— 

cos 9i=£, ipi = 7 t£', (Al) 

where £, £' are some random numbers from the interval 
(—1,1). This algorithm becomes inefficient at low tem¬ 
peratures or strong anisotropy. In this case the majority 
of randomly chosen spin directions has to be rejected due 
to large energy change AH. We use such an update to 
check the results of the second algorithm with tuned step 
of change in the spin direction. 

Our main algorithm for spin change was an algorithm 
where a new spin direction is chosen within a small angle 
near a given spin S 2 65 ! First, a random unit vector w 
perpendicular to the chosen spin S, is generated. Then 
new trial configuration is chosen as 

S' = cos 8i Si + sin 0, ; w, (A2) 

where 9i : the rotation angle from S, to S', is chosen ac¬ 
cording to 

COS 0i — 1 T £(cOS 9 rn ax 1); (A3) 


£ is a random number varying in the interval (0,1), 
the angle 0 max is a maximum allowed amplitude for 
the change of the polar angle Qi of the initial spin S i; 
0 < 9 max ^ 7 r. The new spin direction S' lies within 
a cone around the initial direction with aperture an¬ 
gle 2 9 max and all the directions inside this cone can 
be reached with the same probability®]. The value of 
9 m ax is adjusted after one full MC sweep over the lat¬ 
tice to keep, when possible, the number of accepted spin 
changes around ~ 50%. Also we kept the lower bound 
for 9 max > 7r/6 to prevent too small MC moves which 
are inefficient to thermalize the system. 


This algorithm is not valid at low temperature^^ or 
strong anisotropy, K 1, when the system tends to 

the Ising limit which does not allow the spin flips. We 
improve the situation allowing spins to flip with a certain 
small probability (< 0.1—0.2). With this modification we 
reproduce the correct values of the critical tem perature 
T c for the Heisenberg, T c ~ i .44 HUMIZHSHSi] and the 
Ising, T c ~ 4.51, models! 6 ^ 
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